Dynamical mean-field theory from a quantum chemical perspective 
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We investigate the dynamical mean-field theory (DMFT) from a quantum chemical perspective. 
Dynamical mean-field theory offers a formalism to extend quantum chemical methods for finite 
systems to infinite periodic problems within a local correlation approximation. In addition, quantum 
chemical techniques can be used to construct new ab-initio Hamiltonians and impurity solvers for 
DMFT. Here we explore some ways in which these things may be achieved. First, we present an 
informal overview of dynamical mean-field theory to connect to quantum chemical language. Next we 
describe an implementation of dynamical mean-field theory where we start from an ab-initio Hartree- 
Fock Hamiltonian that avoids double counting issues present in many applications of DMFT. We 
then explore the use of the configuration interaction hierarchy in DMFT as an approximate solver 
for the impurity problem. We also investigate some numerical issues of convergence within DMFT. 
Our studies are carried out in the context of the cubic hydrogen model, a simple but challenging 
test for correlation methods. Finally we finish with some conclusions for future directions. 
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I. INTRODUCTION 

In molecular quantum chemistry, the use of systematic 
hierarchies of electron correlation methods to obtain con- 
vergent solutions of the many-electron Schrodinger equa- 
tion has proven very successful. For example, the hier- 
archy of second-order Moller-Plesset perturbation theory 
(MP2), coupled cluster singles doubles theory (CCSD), 
and coupled cluster singles doubles theory with pertur- 
bative triples (CCSD(T)) can be used (when strong cor- 
relation effects are absent) to obtain properties of many 
small molecules with chemical accuracy 11. The compu- 
tational scalings of the above methods are respectively 
n 5 , n e , and n 7 , where n is the size of the basis, which 
seems to limit them to very small systems. However, lo- 
cal correlation techniques can further be used to reduce 
the above scalings in large systems to n, and this has 
extended the applicability of such quantum chemical hi- 
erarchies to systems with as many as a thousand atoms 

Less progress has been made, however, in the use of 
such quantum chemical hierarchies in infinite systems 
such as crystalline solids. We recall briefly the reasons 
why. Consider a molecular crystal, where the molecular 
unit cell is represented by a basis of n orbitals. Assuming 
V cells in the Brillouin zone of the crystal, the solid is 
then represented by a basis of nV orbitals. In density 
functional theory (computationally a single-electron the- 
ory) the cost of the calculation scales as the third power 
of the number of orbitals. However, translational symme- 
try means that one-electron operators (such as the Kohn- 
Sham Hamiltonian) separate into V blocks along the di- 
agonal, and the crystal calculation can be performed for 
only V times the cost of the molecular calculation, rather 
than V 3 times, if translational symmetry were absent. In 
correlated calculations, translational symmetry yields a 
less dramatic advantage. For example, for second-order 
Moller-Plesset perturbation theory, while the molecular 



calculation scales as n 5 , the scaling of the crystal calcu- 
lation with translational symmetry is n 5 V 3 , and there is 
still a very steep and prohibitive cost dependence on the 
size of the Brillouin zone |6j. 

Locality of correlation suggests that a formal high scal- 
ing with Brillouin zone size can be avoided in physical 
systems. (Indeed there are many current efforts under- 
way to explore local correlation methods in the crystal 
setting) 0, || . We can then imagine starting with a dif- 
ferent picture of a crystal which is more local in nature. 
Consider a unit cell in a crystal. It is embedded in a 
medium, namely, the rest of the crystal. Translational 
symmetry implies that the medium consists of the same 
unit cells as the embedded cell, and thus an appropriate 
embedding theory for a crystal should take on a self- 
consistent nature. If we were to carry out the embedding 
exactly, we should not expect any less cost than the full 
crystal calculation. However, if we make the assumption 
that we will neglect (in some manner) inter-cell correla- 
tions due to locality, then we can expect the high scaling 
with Brillouin zone size to vanish, since the theory takes 
on the form of a self-consistent theory for a single unit 
cell. 

Recently, dynamical mean-field theory (DMFT) has 
been applied with success to strongly correlated crys- 
tal problems, which are typically not well described by 
density functional theory or low-order Green's function 
techniques foj-TToj] . Note that in this paper, we will use 
the term DMFT in a general sense, to mean not only the 
single-site variant but also its cluster and multi-orbital 
extensions [17j . From one perspective, dynamical mean- 
field theory can be viewed as a framework which realises 
the self-consistent embedding with local correlation view 
of a crystal described above. DMFT is formulated in 
the language of Green's functions, and has the form of a 
self-consistent theory for the Green's function of a unit 
cell (which may be a primitive cell, or more generally a 
computational supercell). The local correlation approx- 



imation is expressed by assuming that the self-energy is 
local i.e. inter-cell elements of the self-energy vanish, 
or in momentum space, that the self-energy is momen- 
tum independent. It is important to note that although 
correlation effects are neglected between unit-cells, one- 
electron delocalisation effects between unit cells are in- 
cluded. This, together with the self-consistent nature 
of the embedding distinguishes the physics contained in 
DMFT from that in simpler quantum chemical embed- 
ding formalisms, such as QM/MM theory [Hj|. DMFT 
has some connections in spirit also to density functional 
embedding methods [ijl [2(| , although the use of Green's 
functions avoids the need to approximate a non-explicit 
non-additive kinetic energy functional. 

There are several ways in which DMFT can benefit the 
traditional quantum chemical correlation hierarchy and 
vice versa. First, DMFT provides a framework through 
which quantum chemical methods for finite systems can 
be translated to the infinite crystal through the local cor- 
relation approximation, avoiding the cost of correlated 
Brillouin zone sampling. (This is true even for non-size- 
extensive methods such as configuration interaction, as 
one is treating the correlation only within a unit cell 
and a bath, not the whole crystal simultaneously). The 
natural way to combine quantum chemical wavefunction 
methods with DMFT is through the discrete bath formu- 
lation of DMFT, where we need to determine the Green's 
function of a unit cell coupled to a finite non-interacting 
bath, a so-called impurity problem. Second, quantum 
chemistry provides systematic ways to treat many-body 
correlations in the DMFT framework. These quantum 
chemical solvers are of a different nature to many of the 
currently used DMFT approximations. Finally, quan- 
tum chemical methods and basis sets allow us to define 
the ab-initio Hamiltonian and matrix elements needed 
to carry out DMFT calculations in real systems, while 
avoiding the empirical parametrisation and double count- 
ing corrections that are currently part of the DFT-DMFT 
framework. 

The current work can be viewed as taking first steps 
along some of the lines described above. We aim to do 
several things in this paper. First, we provide an infor- 
mal description of DMFT from an embedding perspec- 
tive. While we do not introduce new ideas in this context, 
we hope this description may be helpful in forming con- 
nections to quantum chemical approximations. Second, 
we explore quantum chemical wavefunction correlation 
methods (more specifically, the configuration interaction 
hierarchy) in the DMFT framework within the discrete 
bath formulation. These wavefunction methods are used 
as approximate solvers for the DMFT impurity problem. 
Third, we define the DMFT Hamiltonian starting from 
ab-initio Hartree-Fock theory for the crystal, avoiding 
any double counting or empirical approximations. (Here 
we point out Rcf. 21] the preprint of which appeared as 
this work was prepared for submission, which also starts 
from HF theory to avoid double counting, though in the 
different context of DMFT as applied to a finite system) . 



Fourth, we explore some of the basic numerics of the 
DMFT framework, such as the fitting and convergence 
of the finite bath approximation, and the convergence of 
the self-consistency. We explore all these questions in 
the context of a simple model system, cubic hydrogen 
crystal. While a simple system, the correlation in cubic 
hydrogen can be tuned from the weak to strong limit as 
a function of the lattice spacing, and at least in certain 
regimes, contains correlation features (such as the three 
peak structure of the density states in the intermediate 
regime) that to date can only be captured within the 
DMFT framework. 

The structure of the paper is as follows. We begin 
in section [TT] with an overview of the DMFT formalism, 
starting with a recap of relevant theory of Green's func- 
tions, then proceeding to a general discussion of DMFT 
self-consistency and embedding, the formulation of the 
impurity problem and the many-body solver, and the def- 
inition of the DMFT Hamiltonian starting from Hartrcc- 
Fock theory to avoid double counting. Section |HT] sum- 
marises our implementation of the DMFT algorithm. 
Section IIVI describes our exploration of several aspects 
of the marriage of DMFT and quantum chemistry meth- 
ods and DMFT numerics in the cubic hydrogen system, 
including the use of the configuration interaction hier- 
archy as a solver, the convergence of the DMFT self- 
consistency, and the convergence of the DMFT calcu- 
lations as a function of the bath size. We present our 
conclusions in section IVl 



II. AN INFORMAL OVERVIEW OF DMFT 

A. Summary of Green's function formalism 

To keep our discussion self-contained and to establish 
notation, we begin by recalling some of the basic results 
from the theory of Green's functions. More detailed ex- 
position of Green's functions can be found, for example, 
in [22j. Given a Hamiltonian H and chemical potential fi, 
at zero-temperature the Green's function G(lu) is defined 
as 



Gij(uj) = (*oh 



1 



uj + fi - {H - E ) + iO 



1 



w + fi + (H - E ) - iO 



o}l*o> 

ai\*o) (1) 



where i,j label the orthogonal one-particle basis, and Wo 
and Eg are the ground-state eigenfunction and eigenvalue 
of H, respectively. G{uj) explicitly determines many of 
the interesting properties of the system. For example 
the single-particle density matrix P, electronic energy E, 
and spectral function (density of states) A(w) are given 
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respectively by 

/oo 
e lu>0+ G(uj)duj (2) 
- oo 

1 r°° 

E = --i e tuj0 +Tr{(h + u)G(Lj)]du) (3) 

2 J-oo 

A(w) = --SG{uj + iO+) (4) 

7T 

In general, a; is a complex variable. Real u> corresponds 
to physical frequencies, and for example, the density of 
states (UJ is defined on the real axis. However, it is often 
more convenient to work away from the real axis. For 
example, expectation values such as Eqs. @, ©, should 
be evaluated on contours away from the real axis to avoid 
singularities in the numerical integration. 

In a crystal, we assume a localized orthogonal one- 
particle basis of dimension n in each unit cell. Us- 
ing translational invariance, it is sufficient to write the 
Green's function as G(R, oj), where R is the translation 
vector between unit cells and for each R, uj, G(R, uj) is an 
n x n matrix. We shall often refer to the Green's function 
of a unit cell in this work as the local Green's function. 
The local Green's function is then the block of G(R, u>) 
at the origin R = and we denote this by G(Ro, uj). The 
local Green's function determines the local observables, 
such as the density matrix of the unit cell, or the local 
density of states, via formulae analogous to Eqs. ©, ((U). 
With periodicity, we can also work in the reciprocal k- 
space. The fc-space Green's function G(k, uj) is defined 
from the Fourier transform 

G(k,w) = ^G(R,w)exp(ik-R) (5) 

R 

and the local Green's function is obtained from the in- 
verse transform as 

G(R ,^) = i^G(k, w ) (6) 

k 

where V is the volume of the Brillouin zone. 

When the finite system Hamiltonian is of single par- 
ticle form, h = hij a t a ji the corresponding non- 
interacting Green's function is obtained from the one- 
electron matrix h as 

g(w) = [(uj + f i + iO ± )l - h}- 1 (7) 

where we use the convention of lower case g(w) and h(uj) 
to denote quantities associated with a non-interacting 
problem, and the infinitesimal broadening 0± is positive 
or negative depending on the sign of uj. In a periodic 
crystal, we obtain the non-interacting Green's function 
in A;-space from the fc-space Hamiltonian h(k) for each k 
point, 

g(k >W ) = [(^ + M + iO±)l-h(k)]- 1 (8) 



Green's functions G(uj), G'(uj) corresponding to differ- 
ent Hamiltonians H, H' are related through frequency de- 
pendent one-particle potentials termed self-energies. The 
self-energy S(w) is defined via the Dyson equation as 

= G /_1 (u;) — G _1 (w) (9) 

It contains all the physical effects associated with the per- 
turbation H' — H. For example, we can exactly relate the 
non-interacting Green's function g(w) from Eq. ([7]) asso- 
ciated with non-interacting Hamiltonian h, and the inter- 
acting Green's function G(uj) associated with interacting 
Hamiltonian H, through a Coulombic self-energy. From 
the explicit form of the non-interacting Green's function 
g(w), the Dyson equation in this case is 

G~ 1 (u;) = (u) + fj, + iO ± )l-h-'E(u)) (10) 

In a periodic system, the above equation holds at each 
k where the self-energy S(k, ui) now also acquires a k- 
dependence, 

G^k, w) = (w + m + *0±)1 - h(k) - S(k, uj), (11) 
and the local Green's function becomes 

G(Ro, «) = - £)[(w + /x + i0±)l - h(k) - S(k, 

(12) 

In general, it is convenient to relax the assumption of 
orthogonality of the one-particle basis, for example, to 
work with an atomic orbital basis. For this, the unit 
matrix 1 in the above formulae should be replaced by a 
general overlap matrix S, e.g. Eq. (TT2l becomes 

G(Ro,«) = y 5}( w + n + iO±)S(k) - h(k) - S(k, 

(13) 

In addition expectation values must be suitably modified. 
For example, the local spectral function A(Ro, uj) is given 
by 

A(R o , W ) = -^3^G(k, w + z0 + )S(k) (14) 

k 

As our calculations in this work use a non-orthogonal 
basis, we will henceforth use expressions with explicit 
overlap dependence. 

B. DMFT equations 

In DMFT, the central quantity is the local Green's 
function G(Ro,w) (the Green's function of the unit cell) 
which is determined in a self-consistent way, including 
the embedding effects of the crystal within a local self- 
energy (correlation) assumption. Here we describe how 
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the DMFT framework and the local self-energy assump- 
tion and self-consistency are established. Of course, we 
recommend that the reader also consult one of the many 
excellent review articles for further discussion and illu- 
mination of the DMFT formalism [|, 

From Eq. (fl"3|) . we observe that G(Ro, oS) can be calcu- 
lated if we have the exact Coulomb self-energy S(k, oS), 
However, determining S(k, uj) requires solving the many- 
body problem for the whole crystal. Thus the idea in 
DMFT is to approximate S(k, ui) by one of its main 
components, the local self-energy H(cj), in essence, a lo- 
cal correlation approximation. Formally, this is the con- 
tribution to the self-energy of skeleton diagrams in the 
Green's function perturbation theory where the Coulomb 
interaction has all local indices, i.e. all indices local to a 
single unit cell. The DMFT approximation neglects the 
fc-dependence of the self-energy. In real-space, this corre- 
sponds to neglecting off-diagonal terms of the self-energy 
between unit cells. The local approximation is plausi- 
ble due to the local nature of correlation, and in fact 
as the physical dimension or local coordination number 
D — >• oo, the approximation becomes exact With 
the DMFT local approximation, the local Green's func- 
tion defined in Eq. ([6]) is simply 



G(R , w) = - X> + M + »0±)S(k) - h(k) - S(«)]- 



(15) 



Now S(w) is formally defined by contributions of only 
the local Coulomb interaction to the local Green's func- 
tion. However, this is still a many-body problem. In 
DMFT, we usually reformulate the determination of 
E(cj) in terms of the many-body solution of an embed- 
ded, or impurity, problem where we view the unit cell as 
an impurity embedded in a bath of the surrounding crys- 
tal. (The impurity nomenclature originates from impu- 
rity problems in condensed matter such as the Kondo and 
Anderson models, which informed some of the early work 
in DMFT). Within this impurity mapping, the many- 
body determination of the Green's function of the em- 
bedded unit cell or impurity Green's function Gi mp (w), 
defines the local self-energy S(w). 

We discuss the impurity problem, and impurity solvers 
to obtain the self-energy, in more detail the next section. 
We focus for now on how the self-consistent embedding 
is established in DMFT. For the theory to be consistent, 
the impurity Green's function (i.e. the Green's function 
of the embedded unit cell in the impurity model) should 
be equivalent to the actual local Green's function of the 
crystal, at least within the local self-energy approxima- 
tion. This means at self-consistency, 



Gi m p(uj) — G(Rq,w) 



(16) 



The embedding to achieve the equality (|16p can be en- 
forced through an embedding self-energy, the hybridiza- 
tion A (a;). The Dyson equation relating the impurity 



Green's function and the self-energy and hybridization is 
then 

G imp (w) _1 = (uj + fi+ iO±)S - h imp - £(w) - A(w) 

(17) 

where hi mp is a one-electron Hamiltonian in the unit cell. 
Once we have solved the many-body impurity problem 
to obtain Gj mp , Eq. p7|) defines the local self-energy 
through 



S(w) = (uj + fi + iO±)S - h rmp - A(w) - G imp (uj 



- 1 
(18) 



The hybridization A(cj) can also be defined through a 
similar equation from the local Green's function, ob- 
tained from Eq. (|15p 



A(w) = (uj 



*0±)S 



E(w) - G(R ,w)- 1 
(19) 



Schematically therefore, for a given hybridization 
A(w), solution of the impurity problem yields Gi mp (ui) 
and the local self-energy S(w) 

. / x impurity solver _ , . , N 

A(ui) 4 G imp (ui) -> E(w) (20) 

while given the local self-energy, Eq. (TT5t yields the local 
Green's function and the hybridization 



E(w) -> G(k, 



G(Ro,w) -> A(w) 



(21) 



Eq. (|2lj) and Eq. (f20f thus form a self-consistent pair 
of equations for the self-energy and hybridization that 
should be iterated to convergence. These are the DMFT 
self-consistent equations. At the solution point, the im- 
purity Green's function and local Green's function, are 
identical as in Eq. (|16l) . 

We note here that the Green's functions 
G(Ro, uj), Gi mp (w), and the self-energy and hy- 
bridisation £(tj), A(w) are smooth functions away from 
the real axis. For this reason, the impurity problem and 
the numerical implementation of self-consistency are 
always considered on the imaginary axis rather than the 
real axis. Once the self-consistency Eq. (|16p has been 
reached on the imaginary axis, analyticity guarantees 
equivalence of the Green's functions in the whole 
complex plane. One can then use the converged A(w) 
(continued to the real axis) to recalculate properties 
along the real axis, such as spectral functions, as needed. 
(Many quantities, such as density matrices, require only 
information along the imaginary axis, however). 

We recap the main physical effects contained within 
the DMFT treatment - local Coulomb interaction effects 
are included in each unit cell and replicated throughout 
the crystal, in a self-consistent way which takes into ac- 
count the embedding of each unit cell in an environment 
of the others. Long-range Coulomb terms are not in- 
cluded in the theory although they can be systematically 
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added. In section III Dl we describe how the long-range 
terms can be treated at the mean-field level. 

Note that we have assumed in the above that we are 
working at a fixed /i. Normally, however, we are inter- 



JV(Rq) = - — / e^°+Tr 

* J — oo 



will change. Thus together with the self-consistency, 
the chemical potential fi must be adjusted such that 
N(R ) = No(R-o). The full DMFT algorithm to do so is 
summarised in section ITO1 

We now turn to consider the many-body impurity 
problem and methods for its solution. 

C. The impurity problem and solver in the discrete 
bath formulation 

The purpose of the impurity formulation is to obtain an 
impurity Green's function Gi mv {u>) and a corresponding 
self-energy that describes the effects of the local 

Coulomb interaction in the presence of the hybridiza- 
tion A(w). In general, due to its many-body nature, 
the impurity problem cannot be solved exactly. The ap- 
proximate method used to solve the impurity problem is 
known as the impurity solver. 

There are two formulations in which an impurity solver 
can work . In the first one the impurity Green's func- 
tion is expressed as a functional integral, and its determi- 
nation is a problem of high-dimensional integration. This 
is typically performed using Monte Carlo methods such 
as Hirsch-Fye 12311 or continuous time quantum Monte 
Carlo methods |24l - [27j . In this formulation, the bath is 
infinite and one does not deal with it explicitly since it 
can be integrated out thus avoiding any bath discretiza- 



The non-interacting bath yields a hybridization A(w) for 
the impurity orbitals of the form 

M") = £— (24) 
p 1 

In general, we assume that physical A(w) can be ap- 
proximately represented in terms of the non-interacting 



ested not in fixed /i, but in some fixed particle number 
of the crystal per unit cell, iVo(Ro). As S(w) changes, 
iV(Ro), the current particle number in the crystal unit 
cell, given by (using Eqs. © and ([15])) 



duj (22) 

I 

tion error. These methods are powerful but suffer in gen- 
eral from a sign problem, as well as difficulties in obtain- 
ing quantities on the real frequency axis (such as the 
spectral function) which requires analytic continuation. 
We will not discuss the Monte Carlo formulations of the 
solver further here, but we refer the reader to an excellent 
review [Hj]. 

The second formulation describes an impurity model 
with an explicit finite, discrete bath. Here the idea is 
to view the hybridization A(cj) as arising from a one- 
electron coupling between the impurity orbitals (orbitals 
of the unit cell) and a fictitious finite non-interacting 
bath. The relevance of the formulation with discrete 
bath here is that the determination of the impurity 
Green's function reduces to the determination of the 
Green's function of a finite problem, and this can be 
tackled using standard quantum chemistry wavefunction 
techniques which avoid the sign problem encountered in 
Monte Carlo based solvers. We can view then such an 
discrete bath formulation as providing a way to extend 
quantum chemical methods for finite systems to treat the 
infinite crystal, within the DMFT approximation of a lo- 
cal self-energy. 

Denoting the local orbitals by i, j, . . ., and bath orbitals 
by p, q . . ., we can write an impurity Hamiltonian for the 
impurity orbitals and the fictitious non-interacting bath 
as 



(23) 

I 

bath by fitting the couplings V and the energies e, and 
this is generally found to be true. This resembles the 
assumption of non-interacting w-representability of the 
density in density functional theory. Fortunately, the 
convergence of (|24|) with respect to the number of bath 
orbitals is quite rapid; one does not need a bath the size 
of the entire crystal to obtain a good representation of 



J 



£ S(k)[(w + n + iO±)S(k) - h(k) - E(w) 

. k 



imp-\-bath 



ij ijkl ip 



Vi P (a\ap + a^flj) 
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the hybridization. (Recall that the fit to the bath is al- 
ways carried out on the imaginary frequency axis, where 
A (a;) is very smooth). 

The form of the bath hybridisation in Eq. (IM1) requires 
that lim^^oo A(w) — > 0. While this is true of physical 
hybridisations in an orthogonal basis, the case of a non- 
orthogonal basis requires a little more care, as discussed 
for example, in Ref. [l5j]. Rearranging Eq. (IT71) and 
inserting the definition of the local Green's function, we 
see that the hybridisation is given by 



A(w) = (u> + n + iO±)Sj 



(25) 



The definition of the impurity overlap matrix Si mp and 
impurity one-electron Hamiltonian hi mp can be viewed as 
adjustable as the equality of the impurity Greens function 
and local crystal Green's function, can be maintained 
through appropriate definitions of the hybridisation and 
self-energy in Eq. (|T7| . Consequently, we choose Si mp 
and hi m p to ensure that the hybridisation can be repre- 
sented by the form Eq. ([24| . Expanding the denominator 
in powers of 1/w, we find that to ensure A(w) vanishes 
like 1/u), we should define the impurity overlap and one- 
electron Hamiltonian as UM 



G - 



(26) 



(27) 



where Eoo = E(oo). 

Now that we have defined a finite Hamiltonian for the 
impurity and a finite bath, the determination of the im- 
purity Green's function Gi mp (uj) is the determination of 
the Green's function of a finite problem. Gi mp {uj) is de- 
fined through Eq. ([T} with the impurity Hamiltonian, 



Gij (u) 
(*o|« 



f 



;4l*o> 



UJ + H — (Himp+bath — Eimp+bath) + i0 3 

<*o|4— -777 " b 1 ^ a *l*o} 



UJ + \l + i^^imp-\-bath ^imp-\-bath ) ^0 



(28) 



where i,j denote the impurity orbitals, i.e. the local or- 
bitals of the unit cell, and Ei mp+ b a th, are the ground- 
state eigenvalue and eigenfunction of Hi mp+ bath- Both 
'to and the corresponding Gj mp (o;) can be determined 
through wavefunction techniques familiar in quantum 
chemistry. 

One subtlety is that the finite problem Vl/o is deter- 
mined for some fixed particle number Ni mp+ b a th (and 
spin, say). In principle, at zero temperature, we should 
use the N^ +bath (and spin) which minimises E imp+bath 
for the given chemical potential This means that we 



have to carry out a search over these quantum num- 
bers. Of course /j, and A(w) are also changing in the 
DMFT iterations, and thus in the discrete bath formula- 
tion, the impurity model is a function of Ni mp+ bath (and 
other quantum numbers), fj,, and A(w). The structure of 
the full self-consistency involving these variables is sum- 
marised in the DMFT algorithm in section Hill 

A popular approach in existing DMFT applications is 
to use full configuration interaction (FCI) called exact 
diagonalization (ED) in solid state physics community 



to solve for ^> and G imp (uj) [11| 



]. From a DMFT 
perspective, the advantage of this approach compared 
to Monte Carlo techniques is that it provides direct ac- 
cess to the calculation of the Green's function on the 
real axis, and consequently the spectral function, with- 
out the need to perform analytic continuation as is used 
in Monte Carlo solvers. In addition, there is no sign 
problem. However, FCI is naturally limited to very small 
numbers of impurity and bath orbitals, and the cost of 
evaluating the Green's function (typically at several hun- 
dred frequencies) means that such calculations are orders 
of magnitude more expensive than typical ground-state 
FCI calculations for molecules. One way to avoid this 
limitation is to employ the various systematic quantum 
chemistry wavefunction hierarchies as impurity solvers. 
We will investigate one such simple approximate solver, 
the configuration interaction hierarchy, in section IIVBI 



D. Eliminating double counting in DMFT through 
Hartree-Fock theory 

In current applications of DMFT to real materials, it 
is common to combine DMFT with a density functional 
derived Hamiltonian, the so-called DFT-DMFT approxi- 
mation [T3, 0, EH • Within this formalism, one does not 
work with a strict ab-initio Hamiltonian, but rather with 
a model Hamiltonian 



H;, 



Hdft 



2 X! w i]k ia\a]aia k - H d . c . (29) 

ijkl£act 



where Hoft is the sum of one-electron Kohn-Sham op- 
erators and Hd.c. is a double-counting correction (see be- 
low). The two electron interaction Wijki is chosen to 
sum over a set of active orbitals in the computational 
unit cell. In transition metal applications, these are 
usually a minimal basis of d or / valence orbitals, the 
idea being that the Coulomb interaction in these orbitals 
should be treated with the explicitly many-body DMFT 
framework, rather than within a DFT functional. While 
Wjiki may be obtained from ab-initio Coulomb integrals 
[30l |31| or derived via e.g. constrained DFT calcula- 
tions iJH, HH, they are best regarded in this approach 
as semi-empirical parameters. The advantage of using 
DMFT in only an active space is that delocalised, itin- 
erant electrons are well treated by existing exchange- 
correlation functionals, and not well-treated within the 
DMFT framework which neglects non-local correlations, 
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while the many-body DMFT framework allows a sys- 
tematic approach to high order strong correlations in 
the d and / shells. The adjustment of Wijki further al- 
lows one to account for effective screening of the active 
space Coulomb matrix elements by long-range correla- 
tions. The DFT-DMFT approach has been successful in 
reproducing many properties of strongly correlated ma- 
terials and an excellent description of the possible appli- 
cations and the way of dealing with the double counting 
correction can be found in Ref. [l2|, EH However, 
there are obvious drawbacks. In particular, the Hamilto- 
nian may be considered to be uncontrolled on two levels. 
Firstly, since exchange-correlation effects in DFT are not 
separated between different orbitals, there is a double 
counting of the Coulomb interaction in Hdft and w. 
This is the origin of the double-counting correction Hd. c . 
which must be adjusted empirically. The double count- 
ing problem is similar to that encountered in molecular 
quantum chemistry when DFT is combined with active 
space wave function methods (35[. Secondly, the use of 
a parametrized form for Wijki must also be regarded as 
unsystematic. 

In the current work, we take a more quantum chemical 
approach to DMFT where we try to retain a strict dia- 
grammatic control over the approximations made. This 
can be achieved by starting with a Hartree-Fock descrip- 
tion of the crystal. Within each unit cell we identify 
an active space, typically a set of localized atomic or- 
bitals. (In fact, in the application to cubic hydrogen 
in this work, all the orbitals in the unit cell will be ac- 
tive). Then, we use DMFT to treat the active space 
Coulomb interaction while the remaining Coulomb inter- 
actions (e.g. long-range Coulomb interactions between 
unit cells, as well the interactions between the active and 
inactive orbitals) are treated through the Hartree-Fock 
mean-field. The Hamiltonian in the active space treated 
within DMFT therefore takes the form 

H imp = (fjj - fij)a\aj + \ Wijkia\a]a k ai 

ijGact ijkl£act 

(30) 

where the /y terms represents the exact subtraction of 
the active-space Hartree-Fock density matrix P HF , con- 
tribution to the mean-field Coulomb treatment 

h= P id F ( w *kij -wukj) (31) 

klGact 

This subtraction exactly eliminates any double count- 
ing between the mean-field and DMFT treatments. Note 
that while the inactive Coulomb interactions (such as the 
long-range Coulomb interactions) are only treated at the 
Hartree-Fock level (which is a severe approximation in 
many solids) the mean-field treatment may be viewed as 
the lowest level of a hierarchy of perturbation treatments 
of these interactions, and is thus systematically improv- 
able. Ref. [2l| . whose preprint appeared as this work 



was prepared for submission, also explores a Hartree- 
Fock starting point to avoid double counting, but in the 
context of DMFT applied to finite systems. 

III. DMFT ALGORITHM 



Algorithm 1 General DMFT loop structure. Note that 
the DMFT self-consistency is carried out on the 
imaginary frequency axis. 
1: for all N imp+bath do 
2: while iV(Ro) / N (R ) do 
3: Choose new \x (e.g. by bisection) 
4: Perform DMFT self-consistency for S(cj), AH- 
5: Calculate E imp+bat h 
6: Calculate JV(R ) 
7: end while 
8: end for 

9: Choose N™p +bath that minimises E imp+bat h 
10: For N™ p+bath and the corresponding fi, A(oj) and im- 
purity model, calculate G(Ro,k>) including quantities on 
the real axis e.g. spectral functions 



Algorithm 2 DMFT self-consistency for A(w), S(w). 
Note that all calculations are done on the imaginary 

frequency axis. 
1: Obtain Hartree-Fock Fock matrix f(k), overlap matrix 

S(k), density matrix P(Ro), and initial guess for A(oj). 
2: while ||E(cj) - £ oid HI > r do 

3: Construct Hamiltonian for impurity orbitals with overlap 

correction (using S(oo)) 
4: Construct bath representation from A(lj) 
5: Calculate impurity Greens function and new self-energy 

E( W ) 

6: Update self-energy E(w), E oid H 
7: Update A(w) 
8: end while 



We now summarize the DMFT algorithm in our cur- 
rent implementation, following the basic ideas outlined 
in the earlier sections. We have implemented our algo- 
rithm in a custom code that interfaces to the Crystal 
Gaussian based periodic code [36[ as well as the Dalton 
molecular code [37| . We recall that within the formula- 
tion with discrete bath, the impurity model is defined as 
a function of three variables: Ni mp+ b at h (particle num- 
ber of the impurity model), /i (chemical potential), and 
the hybridisation A(w) which defines a bath parametri- 
sation. All three have to be determined self-consistently 
together. At the solution point of the DMFT algorithm, 
Nimp+hath minimises the ground-state energy of the im- 
purity model Ei mp+ hath (section III C|) , /i yields the cor- 
rect particle number per unit cell of the crystal iV(Ro) 
(Eq. [22]), and A(u>) satisfies the DMFT self-consistency 
conditions (|20p . (f2"Tj) . The high-level loop structure of the 
algorithm is summarised in algorithm [TJ The individual 
steps are 
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1. Loop over possible particle numbers N imp+ bath 
of the impurity model (to determine Ni mp+ b a th 
which minimises the impurity model energy 
E im p+bath(N im p +b ath))- (In principle we should 
search over spin, but we do not do this is in general 
in our applications here). 

2. For each Ni mp+ bath, search over chemical potential 
fi (e.g. by bisection) to satisfy the crystal unit cell 
particle number constraint iV(Ro) = A^(Ro). 

3.-6. For given [i, N imp+ b a th, carry out the DMFT self- 
consistent loop to determine X(w),A(w) and the 
impurity ground state energy Ei mp+ b a th . Note that 
all calculations are here done on the imaginary fre- 
quency axis. 

9.-10. Determine Ni mp +bath which led to the lowest 
Eimp+bath- Using the corresponding \i and hy- 
bridisation parametrisation, which satisfy the crys- 
tal particle number constraint and the DMFT self- 
energy self-consistency equations, recalculate the 
local Greens function G(R ,w) and other desired 
observables, e.g. the local spectral function A(ui) 
along the real axis. 

The DMFT self-consistent loop for S(w), A(w) con- 
stitutes the core part of the algorithm. It is summarised 
in algorithm^ The individual steps are 

1. Initialisation. Perform a Hartree-Fock (HF) cal- 
culation on the crystal in a local basis. Extract 
the converged Fock matrix f (k) and overlap ma- 
trix S(k) in fc-space, and the Hartree-Fock unit- 
cell density matrix P(Ro). The fc-space Fock and 
overlap matrices are then used to construct their 
real-space analogs in the unit-cell. 

2. Begin DMFT self-consistent loop until convergence 
in the self-energy (to within a threshold r) is 
reached. 

3. Impurity Hamiltonian construction. Construct the 
impurity orbital part of the Hamiltonian. The two- 
body integrals Wijki are computed in the same local 
basis as used in the crystal calculation. The one- 
body Hamiltonian for the impurity orbitals hi mp 
is defined as in Eq. ([31]) using the exact sub- 
traction of the mean-field Coulomb treatment i.e. 
himp = f (Ro) — f (Ro), while the overlap of the im- 
purity orbitals is taken as the overlap in the unit- 
cell, S mp = S(R ). Finally, h imp and S imp are 
corrected as in Eqs. (gU), (|2"7)) . 

4. Bath construction. From the hybridisation A(w), 
obtain the bath Hamiltonian parametrisation by 
fitting. In the first iteration, the hybridisation is 



fitted to the Hartree-Fock hybridisation, defined as 

A hf (uj) =(u + /i + iO±)S imp - h imp + (32) 

I -l 

ij>' + A' + iO±)S(k)-f(k) 



This provides a good guess for the DMFT algo- 
rithm. Further details of the bath fitting algorithm 
are given in section IIVDI and in the appendix. 

5. Calculate the ground-state wavefunction of the im- 
purity problem (for given N imp+bath ). Then cal- 
culate the impurity Green's function on the imag- 
inary axis using a truncated configuration interac- 
tion solver, described in section HVBI 

6.-7. Update the self-energy and hybridisation 

A(w) defined through Eqs. (|9]). ([19]). For bet- 
ter convergence, the self-energy is updated in a 
damped fashion, S(w) (1 - a)S(w) + aS^fw), 
where < a < 1. 



IV. BENCHMARK DMFT STUDIES 

We now proceed to our benchmark DMFT studies. In 
particular, we investigate 

1. the preliminary combination of quantum chemical 
and DMFT ideas, using the configuration inter- 
action (CI) hierarchy as a solver for the DMFT 
impurity problem (or conversely, using DMFT to 
extend truncated CI variants to treat the infinite 
crystal), starting from an ab-initio Hartree-Fock 
DMFT Hamiltonian, 

2. the numerical behaviour of the DMFT algorithm, 
including convergence of the self-consistency cycle, 
fitting the hybridisation by a finite bath, and con- 
vergence of correlated properties (such as spectral 
functions) as a function of bath size. We should 
stress that similar studies were caried out before 
using full configuration interaction (FCI) called ex- 
act diagonalization (ED) in solid state physics com- 
munity. Here, however, we will focus on using the 
truncated version of configuration interaction as a 
solver that was developed by us and examine with 
it the questions of interest concerning the numerics 
of the DMFT algorithm. 

Our studies are carried out on an idealised test sys- 
tem, namely (three-dimensional) cubic hydrogen. Hy- 
drogen clusters in 1, 2, and 3-dimensions have been pop- 
ular models in the study of correlation effects in quan- 
tum chemistry, as the correlation can be tuned from the 
weak to the strong regime as the lattice spacing is in- 
creased [H, HI]. Here we study only cubic hydrogen (i.e. 
three dimensions). We use a minimal basis (STO-3G) 
and a unit cell with a single hydrogen atom, and the 
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initial Hartree-Fock crystal calculations are carried out 
using the Gaussian based periodic code Crystal (36[. 
The use of a Gaussian basis means that we employ the 
general non-orthogonal formulation for the Green's func- 
tion quantities in section Hi Al as well as the overlap cor- 
rections to the impurity model Hamiltonian and overlap 
in section III CI Note that the impurity problem in this 
case has only a single Is impurity orbital, and the local 
Green's function also only has a single orbital index. 

We begin with a brief overview of the properties of 
the DMFT solution of the cubic hydrogen model before 
proceeding to discuss the areas above. 



A. The cubic hydrogen solid model 

We have carried out DMFT calculations on the cubic 
hydrogen model for a variety of lattice constants. We find 
that cubic hydrogen exhibits three electronic regimes as 
a function of lattice spacing which are well-known from 
analogous DMFT studies of Hubbard models [13, HH E3, 
Eol . l4l| . We first summarise the main features of the 
spectral functions and the impurity wavefunctions. (The 
spectral functions plotted here are defined as the trace of 
the local spectral function in Eq. (fT4|) ). The regimes are 

• Metallic regime. This occurs with lattice constants 
near equilibrium, and is illustrated by calculations 
at lattice constant 1.4 A. The spectral function dis- 
plays a single broad peak, indicative of metallic be- 
haviour and the delocalised character of the elec- 
trons (Fig. [T]). The metallic nature is also reflected 
in the ground-state wavefunction of the impurity 
model, which is primarily a single determinant, as 
seen from the natural orbital occupancies (Table l!!) 
and from the impurity wavefunction determinant 
analysis (Table [XJ> . Compared to the restricted HF 
spectral function, the correlated DMFT spectral 
function in Fig. [T] displays additional features at 
large frequencies and is broader, but the spectra are 
similar as expected in the weakly correlated regime. 

• Intermediate regime. At intermediate lattice con- 
stants (e.g. 2.25 A and 2.5 A) the spectral func- 
tion develops a three peak structure with fea- 
tures of both the metallic and insulating regime 
(Fig. [[]). In early DMFT work on the Hubbard 
model the central peak was a correlated feature of 
the spectrum not predicted in mean-field theories 
[Tol [Til |40| . The two outer peaks are shifted from 
the ionisation potential and electron affinity of the 
atom. Analysing the impurity wavefunction, we 
find that at both 2.25 A and 2.5 A lattice constants, 
the wavefunction has multideterminantal character 
with significant mixing of open-shell singlets and 
doubly excited determinants into the ground-state 
(see Tables HUH. 

• Mott insulator regime. This occurs at large lat- 
tice constants when the hydrogen atoms assume 



distinct atomic character. This is illustrated by cal- 
culations at lattice constant 6.0 A. (In this limit, 
the DMFT approximation of a local self-energy be- 
comes exact). The spectral function (Fig. [1} dis- 
plays an insulating gap and peaks centered at the 
electron affinity and ionization potential of the hy- 
drogen atom. The impurity wavefunction is a mix- 
ture of open-shell singlets (see Table H]) . We find 
that the singly occupied impurity natural orbitals 
(Tablein]) are respectively localised on the impurity 
and the bath, thus we characterise the impurity 
ground-state as an impurity-bath singlet. (Note 
that the RHF spectral function stays metallic. An 
unrestricted mean-field calculation would yield two 
peaks similar to the DMFT spectral function, but 
at the expense of breaking spin symmetry). 



B. A configuration interaction impurity solver 



As described in section Hi CI once the impurity model 
Hamiltonian has been defined, we can determine the im- 
purity Green's function within a wavefunction formalism. 
Here we investigate the use of the configuration interac- 
tion (CI) hierarchy to construct impurity solvers. We 
can also see this as using the DMFT framework to ex- 
tend configuration interaction to the infinite system. To 
the best of our knowledge, truncated configuration inter- 
action has not previously been explored in the DMFT 
literature, although full configuration interaction (exact 
diagonalisation) has been widely used [ll|,|42|. By consid- 
ering CI at an arbitrary excitation level we obtain a hier- 
archy of impurity solvers that can, with increasing effort, 
be systematically converged to the exact full CI limit, 
within the given bath parametrisation. We have based 
our implementation on the arbitrary excitation level CI 
program in Dalton 37[ . Our code allows the additional 
possibility of defining restricted active spaces [43( . How- 
ever, for the simple cubic hydrogen model, we find that 
the restricted active space methodology is not necessary. 
Detailed studies of the active space flexibility of the solver 
will thus be presented elsewhere. 

To carry out CI we define a starting determinant in a 
"molecular orbital" basis. Note that this is quite different 
from how exact diagonalisation is used in DMFT, where 
the one-particle basis is chosen to simply be the site basis 
(atomic orbital basis) of the impurity and the bath. Of 
course, the result of exact diagonalisation is independent 
of the choice of one-particle basis, and in model prob- 
lems (such as the Hubbard model), the Hamiltonian has 
a particularly simple local form in the site basis of the 
impurity and bath. However, for truncated configuration 
interaction the choice of starting orbital basis is of course 
much more important. Here we take the molecular or- 
bitals to be the eigenfunctions of the Fock operator of 
the impurity and bath Hamiltonian Himp+bath, Eq. (|2"5]l 
(this is obtained by replacing the impurity part of the 
Hamiltonian by the impurity Fock operator f appearing 
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FIG. 1: Spectral functions (density of states) from FCI, CISD, and RHF calculations for cubic hydrogen, at various lattice 
constants. 

A) ao = 1.40 A, 9 bath orbitals, 300 frequency points. 
C) ao = 2.50 A, 9 bath orbitals, 300 frequency points. 



B) ao = 2.25 A, 9 bath orbitals, 300 frequency points. 
D) ao = 6.00 A, 9 bath orbitals, 300 frequency points. 
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TABLE I: Total weight of CI coefficients of different classes of determinants (Hartree-Fock (HF), singly-excited (S), doubly- 
excited (D)) in the ground-state wavefunction of the impurity model as a function of lattice constant ao. 



excitation level 


a = 1.4 


a = 2.25 


ao = 2.5 


ao = 6.0 


HF 


0.880 


0.755 


0.676 


0.034 


S 


0.040 


0.014 


0.087 


0.941 


D 


0.000 


0.000 


0.098 


0.000 


number of dets with cf > 0.01 


5 


8 


6 


5 




0.920 


0.769 


0.861 


0.975 



in Eq. (|30p ). From the lowest energy orbitals we then 
populate a ground-state determinant and define the set of 
singles, doubles, and higher excited determinant spaces 
as in a conventional CI approximation. We calculate the 
ground-state impurity wavefunction within the given CI 
space, generating a CI vector if> and a ground-state en- 
ergy Eimp+bath- We then evaluate the Green's function 
(128p by solving the two intermediate linear equations for 



X», and Xj 

[(w + jU + Eimp+bath)!- — H imp+ b at / t )]X i (w) = B.j, (33) 

[(w + H — E imp+ bath)l + Himp+bath)]Xj (uj) = Bj, (34) 

B j = Crf 

where Cj, Ci, Himp+bath are representations of the im- 
purity orbital creation, annihilation operators and impu- 
rity and bath Hamiltonian operator in the truncated CI 
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TABLE II: Impurity model natural orbital occupancies for cu- 
bic hydrogen (9 bath orbitals) as a function of lattice constant 
a . 







level 


1-3 


4 


5 


6 


7 


8-10 


a a 


= 1.4 


FCI 


2.000 


1.999 


1.905 


0.095 


0.001 


0.000 






CISD 


2.000 


1.999 


1.905 


0.095 


0.001 


0.000 


a 


= 2.25 


FCI 


2.000 


1.998 


1.718 


0.282 


0.002 


0.000 






CISD 


2.000 


1.999 


1.720 


0.280 


0.001 


0.000 


a 


= 2.5 


FCI 


2.000 


1.999 


1.528 


0.472 


0.001 


0.000 






CISD 


2.000 


1.999 


1.531 


0.469 


0.001 


0.000 


ao 


= 6.0 


FCI 


2.000 


2.000 


1.000 


1.000 


0.000 


0.000 






CISD 


2.000 


2.000 


1.000 


1.000 


0.000 


0.000 



space. (Note, for the N + 1 and N — 1 particle spaces 
accessed by the creation and annihilation operators, we 
consider the space of all determinants that are connected 
to the N particle truncated CI space for the ground-state 
calculation) . u) can be either purely imaginary (as used in 
the DMFT self-consistency cycle) or it can be real, with 
a small imaginary broadening it], when calculating the 
spectral function. The Green's function matrix element 
is then obtained via 

Gij = BiXj + BjX,; (35) 

The solution of the linear equations (|35p can be achieved 
via a variety of iterative algorithms. Our implementation 
follows the algorithm for CI response properties described 
in Ref. [Hj], adapted to truncated CI spaces. 

Our calculations have demonstrated that in the molec- 
ular orbital basis the modest variant of truncated config- 
uration interaction, namly CISD, where the Hilbert space 
is truncated to contain only singly and doubly excited de- 
terminants, was completely sufficient to illustrate all the 
regimes of the hydrogen solid. In Fig. [1] and Table HH we 
show the CISD and FCI local spectral function and im- 
purity natural orbital occupations in the three electronic 
regimes of cubic hydrogen. In the metallic regime, the 
CISD spectral function is completely indistinguishable 
from the FCI spectral function, and the same is true for 
the impurity orbital natural occupation numbers. In the 
intermediate regime, for the lattice spacings 2.25 A and 
2.5 A we expect correlation effects to be stronger. How- 
ever, the impurity natural orbital occupations show that 
there are only two natural orbitals with significant partial 
occupancy, and thus CISD is a very good approximation 
to FCI. This is reflected in both the spectral functions 
in Fig. [T] where CISD and FCI agree very well, as well 
as in the natural orbital occupation numbers, although 
CISD is not as close an approximation in this case to 
FCI as it is in the metallic regime. Finally, in the Mott 
insulator regime, the analysis of the occupation numbers 
shows again that there are only two orbitals with signifi- 
cant partial occupancies and the FCI and CISD spectral 
functions and impurity natural occupation numbers are 
again indistinguishable. 

The near-exactness of the CISD level of impurity solver 
is a feature of the simplicity of the cubic hydrogen model 



system but also reflects the compactness of the CI expan- 
sion when one is using an appropriate one-particle start- 
ing basis, in this case the molecular orbital basis rather 
than the site basis. We expect that more complex solids 
will pose greater challenges and require higher levels of 
excitation in the configuration interaction solver, and 
these issues will be examined elsewhere. Nonetheless, the 
good performance of the single and doubles level trunca- 
tion suggests that it will be promising to explore system- 
atic wavef unction hierarchies in more complex problems, 
which may be infeasible in the exact diagonalisation ap- 
proach. 



C. DMFT numerics: self-consistency 

As discussed in our overview of DMFT and our specifi- 
cation of our implementation in section Hill the impurity 
model particle number Ni mp+ b at h, chemical potential /i, 
and hybridisation A(u>) and self-energy S(w) must all be 
determined self-consistently. The determination of the 
optimal impurity model particle number and chemical 
potential are discrete and continuous searches over single 
variables which are essentially robust. In contrast, the 
self-consistency condition for A(w) and S(w) are multi- 
dimensional equations. Here we examine the convergence 
of the self-consistency cycle for the self-energy S(cj) in 
the loop given by steps 3.-6. in algorithm [2] 

In Fig. [5] we examine the spectral functions obtained 
at the CISD level in the three electronic regimes of cu- 
bic hydrogen as a function of the number of iterations 
of the self-consistency cycle. Generally, we find that con- 
vergence is very rapid. In the case of the metallic regime, 
the spectral function appears to converge after 5 itera- 
tions. In the intermediate regime, for lattice constant 
a = 2.25 A 

the spectral function also converges after 2 iterations. 
At the slightly larger lattice constant ao = 2.5 A, conver- 
gence is a little slower and the spectral function requires 
4 iterations to converge. Finally, as we enter the Mott 
insulating regime, convergence is once again rapid and 
the spectral function converges after 2 iterations. 

The same convergence behaviour is observed in the 
electronic structure of the impurity problem. In Table Hill 
we show the natural orbital occupation numbers of the 
impurity problem corresponding to ao = 2.5 A. These 
numbers were obtained using the CISD solver. (Addi- 
tional tables corresponding to the other lattice constants 
are given in the supplementary material [45]). We see 
that convergence in the 2nd decimal place is reached af- 
ter 5 iterations. 

Overall, we find that at least for the spectral functions 
of the cubic hydrogen model, only a few iterations of 
self-consistency are already sufficient. For quantitative 
properties, such as total energy evaluation of the total 
energy with chemical accuracy, we expect, however, to 
need a tighter convergence. 
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FIG. 2: Spectral function (density of states) obtained with CISD as a solver during the iterations of the self-consistency cycle 
for cubic hydrogen, at various lattice constants. 

A) ao = 1.40 A, 9 bath orbitals, 300 frequency points B) ao = 2.25 A, 9 bath orbitals, 300 frequency points 
C) ao = 2.50 A, 9 bath orbitals, 300 frequency points D) ao = 6.00 A, 9 bath orbitals, 300 frequency points 
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TABLE III: Natural orbital occupancies obtained with CISD 
solver during the iterations of self-consistent cycle for cubic 



hydrogen, ao = 2.5 A, 9 bath orbitals, for exact parame- 
ters used to converge self-consistency see supplementary ma- 
terial ftfj. 

iter/orb no. 1-3 4 5 6 7 8-10 

1 2.000 1.999 1.720 0.280 0.001 0.000 

2 2.000 1.998 1.583 0.417 0.002 0.000 

3 2.000 1.999 1.556 0.444 0.001 0.000 

4 2.000 1.999 1.543 0.457 0.001 0.000 

5 2.000 1.999 1.537 0.463 0.001 0.000 

6 2.000 1.999 1.533 0.467 0.001 0.000 

7 2.000 1.999 1.531 0.469 0.001 0.000 



D. DMFT numerics: convergence with bath size 

As discussed in section Hi CI when dealing with an ex- 
plicit bath the hybridisation A(o>) is parametrised by a 
finite bath, and all quantities must then be converged 
with respect to the number of bath orbitals. There are 
two aspects of bath convergence to explore. How difficult 
is the numerical problem of fitting the hybridization to 
the bath couplings e p and V V {1 How rapidly do the rel- 
evant correlated quantities (such as the DMFT spectral 
functions) converge with bath size? In the latter case, the 
ability of the truncated configuration interaction solver 
(here CISD) introduced in section ITVBI to access larger 
bath sizes than available to exact diagonalisation, pro- 
vides a new capability to examine bath convergence. 

We first discuss the numerical fitting and quality of 
representation of the hybridization A (a;) as a function of 
the number of bath orbitals with couplings e p and V P i. 
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We determine the bath parameters e p and V P i by fitting 
A (a;) to the form (|24|) . In principle, one could carry 
out the fit using any set of frequencies, but following 
standard practice, we fit along the imaginary frequency 
axis, where the hybridisation is a smooth function, and 
use an equally spaced set of frequencies ui n (Matsubara 
frequencies) 



(2n+ 1)tt 
u n = - j^, n = 0,l,2. 



(36) 



where f3, the inverse temperature, determines the spac- 
ing. The choice of j3 is somewhat arbitrary, but to repro- 
duce spectral functions over a given range of frequencies, 
we find that it is reasonable to take /3 to correspond to a 
similar range of frequencies on the imaginary axis. 

Fitting to Eq. (f2~4")l is a highly nonlinear fit. We find 
that the final fit quality depends strongly on the initial 
choice of the parameters. We have established an ini- 
tialisation procedure to obtain a reasonable set of start- 
ing e p and Vpi, described in the appendix. From this 
initial set, we use a Levenberg-Marquadt algorithm to 
minimise the metric Yl n ij l^iji^n) — Ay'(u>n)| to refine 
the bath parameters. As described in section III CI the 
non-orthogonal orbital corrections for the impurity over- 
lap and Hamiltonian (|26l) . (l27|) , are essential for obtain- 
ing a reasonable fit when the underlying crystal basis is 
non-orthogonal. However, we find also that if we arti- 
ficially set the overlap matrix S(k) to the unit matrix, 
and proceed to fit the hybridization functions obtained 
in this way, considerably better fits are easily obtained. 
This suggests that it will be more efficient in the future 
to work within a local orthogonal basis for the crystal, 
rather than the Gaussian basis currently used. 

We show the results of the fitting procedure for the real 
and imaginary parts of the Hartree-Fock hybridization 
(defined in section in the metallic regime in Fig. [3] 
and Fig. [U Similar studies of illustrating difference be- 
tween Green's functions obtained for different number of 
bath orbitals can be found in Appendix C of Ref. [ll[ 
or for cluster DMFT in Ref. [13 • It is evident that 
the fit becomes better as we increase the number of bath 
orbitals, and indeed with 5 bath orbitals the fits appear 
exact to the eye. However, the quality of the fit along the 
imaginary axis does not necessarily guarantee the same 
quality of reproduction of properties along the real axis. 
In Fig. [5]we show the convergence of the accuracy of the 
impurity spectral function, — iQTrGi mp (w) to the corre- 
sponding Hartree-Fock quantity — i3Trg(R , lu). (Note 
that this is not the physical local spectral function, which 
must be defined in a non-orthogonal basis which an ad- 
ditional overlap factor, as in Eq. ((14)) ). For two bath 
orbitals the fit on the imaginary axis is poor and the 
spectral function on the real axis is poorly represented 
as well. Once the number of bath orbitals is increased 
to five orbitals, the error of the fit on the imaginary axis 
becomes quite small and the spectral function becomes 
appropriately improved. However, the rate of the im- 
provement of the spectral function with respect to the 



FIG. 3: Fitting accuracy for the real part of the hybridization 
Re(A.(iu))) for various numbers of bath orbitals. The number 
of frequencies employed was 128 and /3 = 128. 
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FIG. 4: Fitting accuracy for the imaginary part of the hy- 
bridization 7m(A(iw)) for various numbers of bath orbitals. 
The number of frequencies employed was 128 and ft = 128. 
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number of bath orbitals is slower than the improvement 
of the fit on the imaginary axis, as it is much less smooth. 
Note that for each of the spectral functions in Fig. [5] we 
have chosen a different broadening parameter r\ to reflect 
the changing bath orbital spacing. 

We now turn to the convergence of the correlated 
DMFT quantities as a function of bath size. The need 
to examine this convergence is an essential feature of 
working within the discrete bath formulation. In Fig. 
[5] we present the cubic hydrogen local spectral functions 
obtained using the CISD method as a solver at lattice 
constant 2.25 A using 5, 9, and 19 bath orbitals in the 
impurity model, the latter bath size being comfortably 
beyond what can be studied using exact diagonalisation. 
In addition, in Table IIVI we also present the impurity 
natural occupation numbers calculated with CISD solver 
with the different bath sizes as a more quantitative test of 
the bath size convergence. Similar studies of the conver- 
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FIG. 5: Fitting accuracy with different number of bath or- FIG. 6: Spectral function (density of states) obtained with 
bitals for the Hartree-Fock impurity spectral function of cubic CISD solver for different number of bath orbitals for cubic 
hydrogen. The number of frequencies employed was 128 and hydrogen, ao — 2.25 A. 
P = 128. 

" SO l 1 1 1 1 1 1 1 1 




TABLE IV: Impurity natural orbital occupancies obtained 
with CISD solver for cubic hydrogen at lattice constants 2.25 
A, using 5, 9, and 19 bath orbitals. 



5 bath 


1 


2 


3 


4 


5 


6 


a = 2.25 


2.000 


1.999 


1.710 


0.290 


0.001 


0.000 


9 bath 


1-3 


4 


5 


6 


7 


8-10 


a = 2.25 


2.000 


1.999 


1.720 


0.280 


0.001 


0.000 


19 bath 


1-8 


9 


10 


11 


12 


13-20 


a = 2.25 


2.000 


1.999 


1.739 


0.261 


0.001 


0.000 



gence of the occupation numbers with respect of to the 
bath size while usin g ex act diagonalization as a solver 
can be found in Ref. |47l 

We see that the spectral functions are in fact quite 
similar between the different bath sizes. Indeed already 
the very small 5 bath orbital result is remarkably similar 
to the 19 bath orbital result. This must be considered 
a feature of the simplicity of the cubic hydrogen model 
which has only a single orbital in the unit cell. Examin- 
ing the impurity model natural occupation numbers we 
also see that all bath orbital sizes yield very similar nat- 
ural occupancies with only very small differences. This 
is promising for future applications as it seems only a 
relatively small number of bath orbitals is necessary to 
obtain a converged result. 



V. CONCLUSIONS 

In this work we have carried out an initial study of 
dynamical mean-field theory (DMFT) from a quantum 
chemical perspective. DMFT provides a powerful frame- 
work to extend quantum chemical correlation hierarchies 
to infinite problems through a self-consistent embedding 
view of the crystal. The basic approximation is one of 
a local self-energy, which is a kind of local correlation 



approximation. 

We have explored several ways in which quantum 
chemical ideas can be combined with the DMFT frame- 
work. First, we start with a Hartree-Fock based DMFT 
Hamiltonian which avoids the double counting problems 
of the commonly employed DFT-DMFT scheme. Sec- 
ond, we have investigated the truncated configuration 
interaction (CISD) as an impurity solver. The CI hi- 
erarchy avoids the sign problem inherent to Monte Carlo 
solvers in DMFT, and allows a systematically improvable 
approach to the exact solution. Conversely, the DMFT 
framework enables even truncated CI to be extended to 
the infinite crystal. In the simple but challenging cubic 
hydrogen model we find that CI at the singles and dou- 
bles level already reproduces the structure of the den- 
sity of states in the various electronic regimes with near 
perfect accuracy. Finally, we have carried out an inves- 
tigation of some numerical aspects of the DMFT proce- 
dure, including convergence of the self-consistent cycle 
and convergence of properties with respect to the bath 
discretisation. We find that modest bath sizes, easily 
accessible to the CI solver, already produce converged 
results. 

These investigations should be viewed as first steps and 
there are many avenues to develop these ideas. For ex- 
ample, the Hartree-Fock starting point in DMFT treats 
long-range Coulomb interactions at only the mean-field 
level, neglecting long-range screening. Quantum chem- 
ical perturbation techniques may be useful in treating 
these additional interactions and may prove complemen- 
tary to current Green's function treatments of screening 
[H |3l| . Also, there is a wealth of quantum chemical wave- 
function approximations that could be combined with 
the DMFT framework, the most obvious example being 
coupled cluster theory, which should prove advantageous 
over configuration interaction as the number of impurity 
orbitals increases. 

Additionally, the main ideas in this work, in particular, 
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the use of quantum chemical Hamiltonians and solvers, 
are not limited to the single orbital DMFT that we have 
used to study cubic hydrogen. Their combination with 
multi-orbital and cluster versions of DMFT [HI BIIlO] 
should be investigated. Finally, the possibility of us- 
ing DMFT in finite systems, either within the standard 
DMFT formalism [HIHi] or through a true finite DMFT 
formalism [2l| , or the use of DMFT ideas with quantum 
variables other than the Greens function are further in- 
triguing possibilities for the future. 
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VII. APPENDIX: GUESS FOR BATH FITTING 

To generate some initial guess bath parameters e p and 
V p i for the bath fitting, we follow the procedure below. 
Let us specialise to the case of a single impurity orbital 
where we can drop the i index. Then the bath parametri- 
sation (|24|) becomes 



where we have assumed V p is real. Viewing l/(uj n — e p ) 
as the elements of a matrix M np = l/(cj n — e p ), the above 
becomes the matrix equation 

A„ = J2 M npW P (38) 
p 

where A„ = A(oj n ) and W p — V p . We can invert this 
equation to obtain the couplings 

n 

where we understand M _1 to mean the generalised in- 
verse in the singular value decomposition sense. There 
are now only two remaining issues. First, we have to 
choose a set of e p to define the matrix M. Second, given 
arbitrary A„, W p is not necessarily positive definite (and 
thus does not necessarily yield real couplings V p ) . We find 
the latter to be a problem particularly when the overlap 
matrix (due to non-orthogonality) is significantly differ- 
ent from unity, which further suggests (as discussed in 
section [IV Dj) that it will be advantageous to work in an 
orthogonal basis in the future. 

In the first case, we take roots of the Legendre polyno- 
mial of order P/2 where P is the number of bath levels we 
wish to fit and map them respectively from the [—1, 1] in- 
terval (associated with the Legendre roots) to [0, oo] and 
[-co, 0] using the transformation 1 — x/(X(l + x)), where 
A is a scaling factor that is optimised to produce the best 
fit. In the second case, we simply take V p = ^(W p 
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VIII. SUPPLEMENTARY MATERIAL 



TABLE VII: a = 6.0 A 



A. Impurity natural orbital occupancies obtained 

with CISD solver during the iterations of 
self-consistent cycle for cubic hydrogen at various 
lattice constants, 9 bath orbitals. 



TABLE V: do = 1.4 A 

iter/orb no. 1-3 4 5 6 7 8 9-10 

1 2.000 1.999 1.908 0.091 0.001 0.001 0.000 

2 2.000 1.999 1.886 0.113 0.001 0.001 0.000 

3 2.000 1.999 1.906 0.094 0.001 0.000 0.000 

4 2.000 1.999 1.902 0.098 0.001 0.000 0.000 

5 2.000 1.999 1.898 0.101 0.001 0.001 0.000 

6 2.000 1.999 1.905 0.094 0.001 0.001 0.000 

7 2.000 1.999 1.905 0.094 0.001 0.001 0.000 
8-20 2.000 1.999 1.905 0.095 0.001 0.000 0.000 



TABLE VI: a = 2.25 A 

itcr/orb no. 1-3 4 5 6 7 8djT 

1 2.000 1.999 1.801 0.199 0.002 0.000 

2 2.000 1.999 1.742 0.258 0.001 0.000 

3 2.000 1.999 1.725 0.275 0.001 0.000 

4 2.000 1.999 1.720 0.280 0.001 0.000 

5 2.000 1.999 1.720 0.280 0.001 0.000 



iter/orb no. 


1-3 


4 


5 


6 


7 


8-10 


1 


2.000 


1.998 


1.160 0.840 0.002 0.000 


2 


2.000 


2.000 


1.000 


1.000 


0.000 


0.000 


3 


2.000 


2.000 


1.000 


1.000 


0.000 


0.000 


4 


2.000 


2.000 


1.000 


1.000 


0.000 


0.000 



B. Calculation details 

• 5 bath dmft self-consistency using CISD solver 
was carried out for 200 imaginary frequencies and 
= 100, the used damping factor was a = 0.4, 
convergence threshold on self-energy r = 0.005. 

• 9 bath dmft self-consistency using both CISD and 
FCI solvers for all the lattice constants was carried 
out for 200 imaginary frequencies and = 100, 
the used damping factor was a — 0.7, convergence 
threshold on self-energy r = 0.005. 

• 19 bath dmft self-consistency using CISD solver 
was carried out for 200 imaginary frequencies and 
= 100, the used damping factor was a = 0.8, 
convergence threshold on self-energy r = 0.009. 
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